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ABSTRACT 


Without nuclear testing, advanced simulation and experimental facilities, such as 
the National Ignition Facility (NIF), are essential to assuring safety, reliability, and 
effectiveness of the nuclear force; these capabilities are invaluable to the nation’s 
Stockpile Stewardship Program (SSP). A significant information gap exists in the 
hydrodynamic response to nuclear detonations that occur near the earth’s 
surface. Numerical simulation methods were used to evaluate the hydrodynamic 
response of earth-like materials and to develop the energy coupling/partitioning 
curve for low yield nuclear detonations close to the earth’s surface. Using LLNL’s 
supercomputers and GEODYN hydrodynamic code, the properties of stress, 
pressure, and energy were evaluated for twelve simulated 2.5kT detonations; six 
above the surface and six below the surface. The results indicate stronger air 
blasts for detonations above or near the surface and that energy coupling into the 
ground changes rapidly with detonation location over a very small range between 
the above-ground and below-ground interface. This work serves to provide a 
baseline model to evaluate stress, pressure, and energy in relation to nuclear 
yield close to the earth’s surface. The results support a better understanding of 
the physics of near-surface detonations and also assist in planning future 
experimental work at NIF. 
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I. INTRODUCTION 


From July 1945 to September 1992, the United States conducted 1,054 
nuclear tests and executed two nuclear attacks. Considering the tests performed 
by other nuclear weapon states, the number of nuclear tests done by the U.S. 
nearly doubles that number. Since the initial scientific discovery of nuclear 
fission and the early development of nuclear weapon technology, the 
international situation has matured considerably, and there has been an 
increasing concern about the potentially harmful environmental and public health 
effects of nuclear testing, in addition to the devastating results of nuclear war. 
Over the decades, there have been many efforts on the part of the United States 
to reduce the threat of nuclear weapons. These efforts culminated in August of 

1995 when, then President Bill Clinton, announced the Comprehensive Test Ban 
Treaty (CTBT), which would institute a worldwide ban on the testing of nuclear 
weapons of any yield. This treaty was signed by the United States in September 

1996 at the United Nations, though the agreement has not yet been ratified by 
the U.S. Senate and, thus, has not officially entered into force. Nevertheless, the 
U.S. and other nations have complied with the major features of the treaty, and a 
de facto ban on nuclear weapons testing has emerged. 

Given the current de facto ban on nuclear testing, in order to remain 
confident in the safety, reliability, security and deterrent effectiveness of the 
existing nuclear weapons arsenal, while retaining proficiency in nuclear 
technology, methods other than nuclear testing have been developed [1]. The 
Science-based Stockpile Stewardship (SBSS) Program was established to: 
increase the understanding of the current nuclear stockpile, predict and evaluate 
potential problems as these nuclear weapons age, fix or re-manufacture weapon 
components, and maintain proficiency and credibility to support the deterrent 
value of potential nuclear weapon use in the future [2]. According to the U.S. 
DOE/NNSA, “Computer models and advanced experimental capabilities that 
provide accurate predictive capability of weapon performance in the absence of 
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nuclear testing became the main goal of the National Nuclear Security 
Administration’s Stockpile Stewardship Program.”[2] The National Ignition Facility 
(NIF) at Lawrence Livermore National Laboratory serves as a vital tool to the 
nation’s Stockpile Stewardship Program. NIF is currently the largest laser 
system ever built in the world. It has the capability to achieve pressure and 
temperature conditions that are replicated only in nuclear weapons, the sun, and 
stars, albeit on a much smaller scale [1]. 

During the era of nuclear weapons testing, there were many experiments 
conducted that measured the hydrodynamic properties and shock response of 
materials above and below ground. There is a wealth of knowledge documented 
from the nuclear tests that were conducted, but the majority of these nuclear 
tests took place either at altitudes well above ground level or at depths of burial 
well below the ground surface. Consequently, there is little data of the 
hydrodynamic response of materials to a nuclear detonation at or near the 
earth’s surface. Given the current ban on nuclear testing, large-scale testing of 
weapons cannot fill such data gaps; however, LLNL’s National Ignition Facility 
offers a unique capability to generate data in this region of interest that cannot be 
duplicated by any other conventional means. Employing the laser technology of 
NIF, along with advanced numerical computation methods, data can be produced 
that give insight to both the hydrodynamic response of materials, and the shock 
physics that takes place during and after a nuclear detonation at or near the 
earth’s surface. The data generated will serve particular interests in the fields of 
nuclear detonation detection and verification, nuclear forensics, and structural 
survivability. The primary objective of this thesis research is to perform a series 
of simulation calculations to analyze and support the planning of NIF experiments 
designed to characterize the shock response of a prototypical material from near¬ 
surface detonations, and to generate an energy coupling plot that shows the 
air/ground energy coupling and the energy partitioning between air and ground 
that results from nuclear detonations occurring at or near the earth’s surface. 
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II. PROBLEM SETUP 


A. DEFINING THE PROBLEM 

The principal factors that distinguish a nuclear detonation from 
conventional detonations are specific energy and yield. The two nuclear 
weapons used in the attacks conducted by the United States on Hiroshima and 
Nagasaki had yields of 12-15 kilotons of TNT equivalent and 20-22 kilotons of 
TNT equivalent, respectively. In a NIF experiment simulating an actual weapon 
detonation, the experimental setup must be planned to properly scale down the 
physical dimensions and characteristics of the event. In order to properly mimic 
the conditions of a low yield nuclear explosion, the parameters of energy, length, 
and time must be scaled down to the lab frame. The National Ignition Facility 
has 192 laser beams each with the capability of delivering approximately 9.4kJ of 
energy for a combined total of 1.8MJ of light energy. NIF can generate a peak 
power of 500 trillion watts, which is 1,000 times the electric generating power of 
the entire country [3]. High-precision lasers, such as these, allow for the 
achievement of specific energies that are near those of nuclear devices. This 
opens up a window of opportunity to investigate the physics of near-surface 
nuclear detonations. For the experiments being considered, only four of the 192 
laser beams will be used, with each contributing only 2.5kJ (i.e., operating at 
roughly 25% capacity), to achieve a total deposited energy of lOkJ focused into a 
2mg cylindrical hohlraum,^ which yields a specific energy of 5MJ/g. This 
cylindrical hohlraum has a height and diameter of about 2mm and a volume 100 
times smaller than that of a pencil eraser. 

B. SCALING 

To scale the experiment down to the laboratory frame, the Glasstone/ 
Dolan cube root scaling law was used: 

^ A small gold cylinder named after the German word for “hollow room” that serves as the 
target for a NIF experiment. 
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Theoretically, a given pressure will occur at a distance from an 
explosion that is proportional to the cube root of the energy yield. 
Full-scale tests have shown this relationship between distance and 
energy yield to hold for yields up to (and including) the megaton 
range. [4] 

For the anticipated experiment relating a 2.5kT full-scale test to a NIF 
energy deposition of 10kJ, a scaling factor of 1000 can be derived as follows: 

{1.5kTI\QU)'"^=[1.5kT*{A.\Ux\(fUlkT)l\QU)f'^^~\mQ 
wherelA:r = 4.184*10‘V 

To elaborate, in an actual 2.5 kT yield detonation, the energy release is: 
2.5kT = 10.46 X 10^^ J = 10.46 x 10^ kJ 

This experimental NIF energy deposition of lOkJ in comparison with the 
detonation release of 2.5kT represents an energy reduction factor of: 

(2.5kT/1 OkJ) = (10.46 x 10^ kJ/10 kJ) ~ 10® 

Applying Glasstone and Dolan’s cube root scaling factor rule, the resulting 
scaling factor is 10^, or a factor of 1000 [4]. Another way to characterize the 
cube root scaling factor rule as applied to distance scaling is to recognize that: 

Da= Dr * (Wa/ Wr)'^^ 

where: 

Da=actual distance 
Dr=reference distance 
Wa=actual yield 
Wr=reference yield 

Therefore, in the case being considered where the energy reduction factor 
is 10®, the distance scaling is a factor of one thousand. A similar calculation is 
also applied to the time parameter. As a result, with both time and distance 
properly scaled, the physical properties such as particle velocity, pressure, 
density, and temperature are accurately reflected without scale factors. The 
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primary properties of interest in this thesis are those of pressure and energy, 
which will be used to construct the energy-coupling curve as the height of burst/ 
depth of burial is varied for a 2.5kT simulated explosion. 

For example, Figure 1 presents the results of a simulation of a 10kJ NIF 
experimental explosion at a height of burst of 3cm (corresponding to a 2.5kT 
nuclear detonation at a height of burst of 30m), in which the pressure wave 
reaching a maximum of .8MPa (8 bar), arrives in approximately 80 microseconds 
at a sensor location (Gauge #2) 15 cm from the explosion. This would scale 
back to the real world reference frame as a shock with the same peak pressure 
of 8 bar arriving in approximately 80 milliseconds at a distance of 150 meters 
away from an explosion of 2.5kT yield at a height of burst of 30 meters. 
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III. PHYSICAL CONFIGURATION 

A. NIF CONFIGURATION 

For a series of anticipated NIF experiments to investigate near-surface 
detonation coupling, four of the 192 laser beams will be used. Each beam 
contributes 2.5 kJ of energy over two nanoseconds that will be deposited in a 
2mg cylindrical hohlraum. The four lasers will combine inside the near-vacuum 
target chamber, as shown in Figures 2 and 3, and focus directly into the gold 
hohlraum target. The gold will essentially vaporize and expand as a gold gas pill. 
This is similar to physical conditions in the detonation of a full scale nuclear 
device where all components are vaporized and contribute to the ensuing fireball. 
The gold hohlraum will be located 30mm above the surface of a cylindrical core 
of a ceramic material called Macor^, which for the purpose of these experiments, 
simulates the solid material of the earth’s surface. The properties of the Macor 
ceramic, some of which are listed Table 1, are similar to the properties of many 
solid earth materials, in particular limestone. The Macor ceramic in the 
experiment serves as a model for the ground response to a nuclear detonation 
over limestone. 


Density 

2520 kg/m^ 

Young’s Modulus 

64.1 GPa 

Shear Modulus 

25.4 GPa 

Poisson’s Ratio 

0.26 

Quasi-static compressive strength 

345 MPa 


Table 1. Physical Properties of Macor glass ceramic (From W.Chen, 1997, 

p.1310) 


2 Macor is a machinable white glass ceramic. 
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Figure 3. Target Chamber Exterior 


The experimental target consists of a hohlraum target (representing the 
scaled down nuclear device) and the pre-fabricated Macor cylinder which will 
include three pressure sensors embedded into the structure in addition to the 
hohlraum, and two additional pressure sensors mounted outside of the Macor 
object in order to measure pressure in the region representing the near-surface 
atmosphere. From the focal point of the laser beam, the pressure gauges^ will 
be located at the straight-line distances shown in Table 2. 


2 PG1 and PG2 represent the two atmospheric pressure sensors external to the Macor object 
while GPG1 through 3 represent the pressure sensors embedded within the Macor object. 
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Pressure Gauge 1 (PG1) 

10cm 

PG2 

15cm 

Ground Pressure Gauge 1 (GPG1) 

5cm 

GPG2 

7cm 

GPG3 

10cm 


Table 2. Distances from target to pressure sensors 


' ‘ Pressure Gauges 1 arxJ 2 


4 converging 
NIF beams 


PG1 PG2 


Ground Pressure Gauges 1. 2. and 3 


Macof 


GPG1 


l _ 

GPG3 

GPG2 ^ 

Figure 4. 

Set Up Drawing After [6] 


Figure 4 shows a diagram of the experiment, self contained in a cylindrical 
casing that will be inserted on the diagnostic instrument manipulator (DIM) 
depicted in Figure 5. The subsurface sensor, GPG2, has a quartz piezoelectric 
element with a dynamic range of 0-2 GPa and a bandwidth of 1-100 MFIz. The 
subsurface sensors GPG1 and GPG3, located at 45-degree angles from the 
laser vector in its focal point, have polyvinylidine fluoride (PVDF) piezoelectric 
elements, which provide longer recording times, and each has a dynamic range 
of 0-2 GPa and a bandwidth of 1-100 MFIz. PVDF is a piezoelectric plastic 
polymer that is commonly used as a shock sensor [7]. The above-surface 
sensors PG1 and PG2 have PVDF elements and have a dynamic range of 0-6 
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MPa and a bandwidth of 30kHz to 3MHz. The Macor cylinder will be inserted into 
the NIF target chamber attached to the diagnostic instrument manipulator (DIM), 
shown in Figures 5 and 6. 



Figure 5. DIM and Chamber 



Figure 6. Close Up View of the DIM 


B. CONFIGURATION OF NIF EXPERIMENTAL DESIGN AND 
CORRESPONDING GEODYN SIMULATIONS 

In order to properly model the experiment simulating a nuclear detonation 
over or in limestone, a combination of materials with similar properties to that of 
limestone and air must be used. With Macor as the solid material, the gas mix 
was selected by empirical methods to mimic as close as possible atmospheric 


10 



density. The required gas mix was calculated to consist of 65% Ne, 20% Ar, 10% 
Kr, and 5% Xe. This mixture yields similar results in terms of transmission and 
absorption of the soft x-rays emitted in order to produce a scaled deposition of 
the x-rays and the “double peak” intensity plot that is characteristic of nuclear 
detonations. The first peak in light intensity is due to the initial reaction and 
expansion of hot gas in the form of a fireball. Quickly, the shock wave from the 
explosion expands past the fireball and the air outside the fireball is compressed 
by the shock wave to the point of becoming opaque to light. As the shock wave 
expands spherically, the pressure front decreases and the light from the fireball is 
once again visible [8]. Table 3 presents a summary comparison of the properties 
of the Macor/artificial atmosphere, and limestone/normal atmosphere which 
demonstrates their similarities. 


PROPERTY 

GAS MIX OVER MACOR 

AIR OVER LIMESTONE 

‘Atmospheric’ density 

.00153 g/cc 

.001204 g/cc 

‘Ground’ density 

2.52 g/cc 

2.59 g/cc 

‘Ground’ sound speed 

4.58 km/s 

3.67 km/s 


Table 3. MACOR and Limestone Comparison After [6] 


To achieve the objective of this thesis, a series of simulations using 
LLNL’s GEODYN4 computer code was carried out to analyze the experimental 
configuration for planning purposes and also to develop response projections to 
establish the ground-coupling curve. The first set of simulations were carried out 
to represent the hydrodynamic response to detonations (and NIF experiments 
representing) above-ground detonations. The second set treated below-ground 
detonations. 

1. Above Ground Configuration 

Figure 7 presents an experimental plan and the corresponding simulation 
layout for the above-ground case. 

4 More information about GEODYN can be found at https://e-reports- 
ext.llni.gov/pdf/326532.pdf. 
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Figure 7. 3cm Height of Burst Calculation Set Up; Z-and R- axis in 

millimeters 

In this portion of the simulations, a simplified 2D axisymmetric geometrical 
representation of the physical configuration was constructed. Figure 7 shows the 
output of LLNL’s graphics software, VisIT^, at time zero, as the initial set up in the 
GEODYN computer modeling code. In the above-ground model setup, there are 
five materials—Macor, the atmospheric gas-mix, gold, “air,”6 and steel. The 
purposes of the first three materials are obvious, but the latter two require 


5 VisIT is an interactive paraliei visuaiization and graphical analysis tool for scientific data 
capable of handling data on the terabyte scale. 

5 The “air” is given an extremely low density to simulate the vacuum conditions, through 
which the lasers propagate, inside the NIF target chamber. 





12 










additional explanation. Within the computer model, there must be a replication of 
the vacuum barrier through which the lasers propagate to deposit the energy 
inside of the gas mix for an above-ground calculation. In this simulation, a hollow 
steel cone is used to represent the vacuum barrier and to keep that vacuum 
separate from the gas mix. (See Appendix A) Essentially, the simulation begins 
after all the energy from the laser is deposited into the gold hohlraum. This was 
achieved by creating a steel cone with an apex at the laser focal point. In the 
simulation, the steel cone is filled with “air,” but given a density of almost zero, so 
that it does not contribute to the forming shock wave. The actual pressure 
anticipated inside of the NIF target chamber is close to 10'® millitorr at the time of 
a target shot. By changing the location of the apex of the steel cone (and the 
attached hohlraum target), the desired height of burst for the calculation is 
achieved. In order to record the changes in the hydrodynamic properties as the 
experiment evolves, marker files were embedded throughout the configuration 
with the primary markers being geometrically located at the same places that the 
subsurface and aboveground pressure sensors are located. Each of the five 
materials is indexed for identification with label values from 0 to 4 and given a 
density and a strength parameter and defined by a series of coordinate points. 
The input file for the 10mm height of burst simulation can be found in Appendix A 
for further clarification. 

2. Below Ground GEODYN Configuration 

The depth of burial simulation set up is similar to that of the height of burst 
simulations, with the exception that there is no need to incorporate the steel 
cone, as the energy is deposited directly into a hohlraum embedded in the Macor 
block. Figure 8 shows the initial set up of the simulation, in VisIT, for the below¬ 
ground case. The subsurface calculations are done from an idealized point of 
view. The Macor is uniformly manufactured and the Macor representation is 
uniform. This is a limitation in terms of representing real-world earthen 
structures, which are far from uniform. The possible addition of cracks, porosity, 

or other heterogeneities within the simulated Macor is not addressed in this work 
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but may be a topic for future study. Additionally, in terms of mimicking a setup 
for a NIF experiment there must be a conduit through which the energy, from the 
laser, is deposited. Essentially this would involve constructing an inverted hollow 
cone with an apex beneath the surface. The amount of energy that remains in 
the region of the cone as opposed to being contained completely in the “ground” 
must also be taken into consideration when calculating the energy coupling. It 
must be stressed that the depth of burial simulations conducted in this work were 
done from an ideal standpoint. This would be the real world equivalent of digging 
a cavity into the ground, placing a 10kJ nuclear device in the cavity, and then 
covering the cavity and packing it to the ground’s original density. These depths 
of burial calculations serve as a baseline used to construct the energy coupling 
curve. 


14 



DB: Header 
Cycle: 0 Time:0 


Pseudocdof 
Van color 

■2.tX» 


- 1.500 


If- 1.000 
-0.5000 


0.000 

Max: 2.000 
Min: 0.000 


Z-Axis 



Figure 8. 


-150 -100 -50 g 50 100 150 

R-Axis 

3cm Depth of Burial Calculation Set Up; Z-and R-axis in millimeters 


15 




THIS PAGE INTENTIONALLY LEFT BLANK 


16 



IV. PHYSICS BEHIND THE CODE 


A. GEODYN BASICS 

In order for a hydrodynamic simulation to have validity, it must be 
constrained so as to obey the laws of physics. The hydrodynamic code used in 
this thesis, called GEODYN, was developed at LLNL and incorporates physical 
models to fully describe a broad range of phenomena including shock and 
thermodynamic behavior. It is an Eulerian code with adaptive mesh refinement 
(AMR). There are two basic types of dynamic codes, Lagrangian and Eulerian. 
A Lagrangian code’s mesh, or background, moves with the material, so no mass 
flows between cells. In an Eulerian code, the mesh or background is stationary 
and the material is allowed to move though stationary cells. An Eulerian code is 
analogous to looking at a dust storm through a window whereas a Lagrangian 
code is analogous to floating along with an individual dust particle in the storm. 
The adaptive mesh means that the code has the ability to vary the level of detail 
of the background. In keeping with the dust storm analogy, the adaptive mesh is 
analogous to having the ability to simultaneously have small microscopic 
windows to capture the details of small grains and larger windows for use when 
detail is not needed. An Eulerian code with adaptive mesh refinement such as 
GEODYN allows for rigorous high numerical resolution in areas in one part of a 
problem and less refined, less sensitive areas in another. 

When dealing with shock physics, high temperatures, and nuclear 
interactions, the laws of physics that constrain the results are numerous. The 
two primary areas of physics that are of concern are thermodynamics and 
continuum mechanics. The GEODYN code simplifies the matter without 
deviation from real experimental results by implementing the various physics 
constraints into a formulation based on both the 1st and 2nd laws of 
thermodynamics, and the laws of conservation of mass, momentum and energy 
[9]. 
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GEODYN is a complex modeling code. In order to provide some 
background information on key physics processes taking place in GEODYN, the 
remaining sections in this chapter, (Sections B and C), provide some 
mathematical derivations of some of the physics involved in terms of 
thermodynamics and continuum mechanics respectively. 

B. THERMODYNAMICS 

A major portion of GEODYN, summarized in Table 4, are the constitutive 
equations that express specific entropy, specific internal energy, stress, and heat 
flux in terms of: temperature (0), change in temperature (A0), deformation 
gradient (F) and the change in deformation gradient (dF). 


Specific entropy (s) 

s=s(0,A0,F,dF) 

Specific internal energy (u) 

u=u(0,A0,F,dF) 

Stress (T) 

T=T(0,A0,F,dF) 

Fleat flux (q) 

q=q(0,A0,F,dF) 


Table 4. Constitutive Equations 


These constitutive equations can be derived from the Clausius-Duhem 
inequality^, which has been derived below for clarification. 

Consider the following five parameters of interest: 

1. Specific entropy s 

2. Mass density of the body p 

3. Internal heat supply/ (mass-time) r 

4. Outward heat flux vector q 

5. Temperature 0 

^ The Clausius-Duhem inequality is used in continuum mechanics to express the second law 
of thermodynamics. 
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The entropy (due to heat transfer into a system) is given by the following 
formula^; 


— = j prl9dv-\- j -qlOda 

Volume Surface 

where dv and da represent the defferentials of volume and area respectively. 

However, when dealing with real, irreversible processes the 
thermodynamic equations become inequalities. So, the change, with respect to 
time of the entropy, would be 



Volume Volume Surface 


Using the divergence theorem (integral over area is equal to the integral of 
the divergence over a volume), the expression becomes: 

— j spdv> j rplOdv- j V»{ql6)dv. 

Volume Volume Volume 


This is the integral form of the Clausius-Duhem inequality. The Clausius- 
Duhem inequality in its most recognized form is in the differential form: 

ds 

— >rie-\lp*V*{qie). 
dt 


This formula states that the rate of entropy change, with respect to time, in a 
system is greater than the rate of entropy change due to heat transfer into the 
system. 


C. CONTINUUM MECHANICS 

The next major portion of GEODYN deals with how a material deforms, or 
flows over time, due to energy input; all while conserving mass, energy, and 


^ The mathematical derivations that follow are summaries based upon detailed presentations 
in [9], 
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momentum. The equations used are often in the form of the convective 
derivative (deformation/ flow equation); for example 

^(#) = 4(#)+V*V(#). 

Dt St 

This equation states that the total change, with respect to time, is equal to 
the change in the fixed frame, plus the change in the moving frame. In effect, 
this is a continuity equation. The parameter denoted by the # sign represents a 
history-dependant scalar parameter such as porosity, temperature, or plastic 
strain and V (capitalized) represents the velocity of the fluid medium. The 
derivation of the generic continuity equation is shown below for the purpose of 
providing clarification for what is taking place mathematically in the GEODYN 
code. 

Let L be some intensive (mass independent) property defined over a 
control volume Q. The parameter Q represents sources and sinks in the fluid. 
The rate of change of L is given by the formula below: 

— j Ldv = - LVnda - \Qdv 

^ Q Surface Q 

where dv and da represent the defferentials of volume and area respectively 


Applying the divergence theorem to the 1st term on the right side of the 
equal sign gives: 


A 

dt 


J LJv = -J V • {LV)dv - J Qdv 

Cl Q. 


Combining the terms on side of the equation gives: 

f(^+v.(Ly)+0dv = o 

The above condition must be true for all dv; therefore, the following 
generic continuity equation results: 
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—+v.(Ly)+G = o 

at 


This equation states that the amount of a conserved quantity, inside a 
defined region, can only change by the amount that passes through the boundary 
of the region in addition to any internal sources or sinks. The conservation of 
mass, energy, and momentum all follow from the above generic continuity 
equation. 

In GEODYN, solid is modeled as an extension of a fluid, with the added 
appropriate strength parameters. A LLNL Mie-Gruneisen equation of stated was 
used to calculate the pressure in the Macor for this simulation. For the gold 
source, a tabular equation of state provided by LLNL was used. 


9 The Mie-Gruneisen equation of state is often used to measure the pressures in shock 
compressed soiids. 
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V. METHOD 


A series of simulation calculations were performed in which the height of 
burst and depth of burial ranged from 5mm to 50mm. A simulation was run at 
5mm, 10mm, 20mm, 30mm, 40mm, and 50mm for both above-ground and 
below-ground cases for a total of 12 simulations. On average, each simulation 
required a total 24 hours of computation, broken into three 8-hour periods on 
LLNL supercomputers. Each computation was considered complete once the 
shock front passed the location of the most distant shock sensor, which on 
average tended to be about 100 microseconds (real time), or about 24 hours of 
computation time. To achieve 12 complete data sets ranging from a depth of 
burial of 5cm to a height of burst of 5cm required a time period of three months. 
After each successful calculation, the data were analyzed using both MATLAB^o 
and LLNL’s visualization software VisIJiT GEODYN has the ability to track 
various parameters during the experiment. The parameters of interest are: 
pressure, density, material type, particle velocity in the r and z coordinate 
directions, sound speed, internal energy, total energy, and temperature. As 
stated earlier, in order to construct data showing the energy coupling between 
ground burst and air shock, the primary parameters of interest are pressure, 
particle velocity, time and total energy. Figures 9 through 13 present the 
pressure wave results for the 10mm height of burst calculation at the five sensor 
positions. 


"'0 MATLAB is a scientific software appiication used for matrix manipulation and plotting that 
can interface with other programming languages such as C, C++, and FORTRAN. See 
http://www.mathworks.com/products/matlab/. 

FOR ADDITIONAL DESCRIPTION OF VISIT, See 
https://wci.llnl.gov/codes/visit/about.html. 
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Height of Burst Plots 



Figure 9. Pressure Wave at Air Gauge 1 for 1cm HOB 



Figure 10. Pressure Wave at Air Gauge 2 for 1 cm HOB 
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Pr*»tu<« at Ground Pr*«$uit Gaug* 1 1cm HOB 



Figure 11. Pressure Wave at Ground Gauge 1 for 1 cm HOB 


Pressure at Ground Pressure Gauge 2 1cm HOB 
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Figure 12. Pressure Wave at Ground Gauge 2 for 1 cm HOB 
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Pressure al Ground Pressure Gauge 3 tcm HOB 



Figure 13. Pressure Wave at Ground Gauge 3 for 1 cm HOB 

The most useful information gathered from these plots is the peak 
pressure and the time of arrival of the peak pressure wave. This information 
served as a check for the calibration of the PVDF shock sensors. For Figures 10 
through 13, after the initial pressure waves any other peaks are the results of 
boundary reflections. In the case of Figure 10 for PG2 a second pressure peak 
is visible at 120 microseconds. This second peak is due to the reflection of the 
initial blast wave from the Macor/gas-mix boundary. For Figures 11 though 13 
the many reflection stress waves are the result of free surface reflections from 
the Macor/gas-mix boundary. Some of the interesting phenomena shown in the 
plots above were made clear when the data were analyzed using LLNL’s 
visualization software, VisIT. For example, in the ground sensors, the stress 
waves in tension were easily identified by stepping through the calculation frame- 
by-frame in time. In the simulation, the Macor is surrounded by a free surface 
and thus tension waves occur at the intersection between the Macor block and 
the gas mix. Figure 14 shows a movie frame of the detonation at 11 
microseconds, with pressure on the left and material on the right, where the initial 
stress wave propagating through the Macor object is visible. Figure 15 shows a 
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movie frame, at 24 microseconds, where the stress waves have reached the 
edges of the Macor object and reflected tension waves are clearly visible. 
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Figure 14. Stress wave propagating through the Macor object 
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Figure 15. Tension waves propagating through the Macor object 
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Since, in the actual experiment, the Macor test object will be held on the 
DIM, the real fluctuation in pressure can be expected to be substantially less. 
The recovery period of the pressure sensor will also eliminate much of the 
reflected response, as the sensors will not be able to record continuously as in 
the simulation. 

Depth of Burial Plots 

The pressure waves for the five shock sensors, for a 30mm depth of burial 
simulation, are shown in Figures 16 through 20. Figure 16 shows that there was 
no air shock wave measured by PG1 for the 30mm depth of burial simulation. 
Figure 17 shows a strange wave profile in the shape of a ramp wave, indicative 
of a supported shock, (which should not be the case here as the lOkJ is 
deposited into the gold hohlraum within a few nanoseconds), for the air pressure 
gauge (PG2). 
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Figure 16. Pressure Wave at Air Gauge 1 for 3cm DOB 
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Prt$tu>« at Piataur* Gaug* 2 3cm OOB 
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Figure 17. Strange Pressure Wave at Air Gauge 2 for 3cm DOB 

Upon initial review of the data, it was unclear how the pressure profile in 
Figure 17 translated to what was happening mathematically in the simulation. 
The wave profile appears to indicate a pressure wave reflection that came from 
the boundary. To verify that this was indeed the case, further calculation was 
conducted in which an additional air sensor (PG3) was set up with an r-z 
coordinate of [25,140], which is closer to the boundary of the setup, to see if the 
pressure peak time of arrival was shorter than that in PG2. The results, shown in 
Figure 18, indicate the same anomaly in the pressure wave at the nominal 3rd air 
pressure gauge (PG3). The peak pressure time of arrival for the nominal PG3 
sensor, arriving at approximately 40 microseconds, was shorter than the same 
anomalous wave in Figure 17 for PG2. This shows that the data is consistent 
with a wave that reflected from the boundary. 
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Pressure at Pressure Gauge 3 3cm DOB 



Figure 18. Strange Pressure Wave at Nominal Air Pressure Gauge 3 

In order to eliminate the reflected wave noise from the depth of burial 
calculations the set up of the problem had to be modified. The boundary 
conditions of the problem are set up so that the system conserves energy, mass, 
etc. This issue did not arise in the height of burst calculation because there is a 
single material, “air,” at the boundary. The issue arose in the depth of burial 
calculations because the boundary has two different materials, both Macor and 
gas-mix. To correct the setup, a thin layer of gas-mix was added around the 
Macor object in the input file and the simulations were run again (See Appendix 
B). Figure 19 shows the modified depth of burial setup. When performing the 
pressure analysis on the modified depth of burial calculations, the noise in PG2 
vanished. Additionally the peak pressures in the three “ground” sensors in the 
Macor object were identical, which shows that the modification of the setup did 
not adversely affect the peak pressure readings. Figures 20 and 21 show the 
pressure wave profiles for PG2 and nominal PG3, which show that little to no 
shock was transmitted from the ground to the gas-mix (the pressure signal is 
down an entire order of magnitude from that in Figures 17 and 18), as expected. 
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Figure 19. Modified Depth of Burial Set Up 


)(lo'® Pressure at Pressure Gauge 2 3cm DOB 



Figure 20. Pressure at PG2 with Modified Set up (no noise from boundary) 
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y in'* Pressure at Pressure Gauge 3 3cm DOB 



Figure 21. Pressure profile for nominal pressure gauge 3 (PG3) 


Pressure at Ground Pressure Gauge 1 3cm DOB 



Figure 22. Pressure Wave at Ground Gauge 1 for 3cm DOB 
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Pressure at Ground Pressure Gauge 2 3cm DOB 



Figure 23. Pressure Wave at Ground Gauge 2 for 3cm DOB 


Pressure at Ground Pressure Gauge 3 3cm DOB 



Figure 24. Pressure Wave at Ground Gauge 3 for 3cm DOB 

As one would intuitively guess, the peak pressures calculated in the Macor 
object are much higher for the depth of burial shots. Again for Figures 22 through 
24, the many reflection stress waves are the result of free surface reflections 
from the Macor/gas-mix boundary. 
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VI. ENERGY COUPLING 


In order to illustrate energy coupling and partitioning, two methods were 
used. The first method demonstrates the coupling by looking at the peak 
pressures for two pressure sensors for varied heights of burst and depths of 
burial. Figure 25 shows three pressure plots for three different heights of burst/ 
depth of burial from air pressure gauge 1 (PG1) superimposed on one another. 
Figure 26 shows three stress plots for three different heights of burst/ depth of 
burial from a ground sensor, (GPG3), superimposed on one another. These two 
sensors, PG1 and GPG3, were chosen because both are located at a distance of 
10cm from the point of detonation. The r-z coordinates for PG1 are (8.6cm, 5cm) 
and the r-z coordinates for GPG3 are (7.07cm, 7.07cm) both yielding a range of 
10cm. The three pressure curves correspond to two heights of burst and one 
depth of burial. The plot in Figure 25 shows that the air blast is strongest for 
detonations that occur near the surface. This may seem counterintuitive. The 
reason that the surface burst yields a stronger air blast is that the ground reflects 
much of the pressure wave, and the shock sensors detect the combined 
pressure of both the initial and the reflected shock. 



tjm« (micro $»cofidi) 


Figure 25. Air blast is strongest for surface burst 
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As reasonably expected, the results of the ground pressure/stress plots 
show that the ground shock is strongest for a buried burst. 



Figure 26. Ground shock is the strongest for a buried burst 

The second method used to illustrate energy coupling and partitioning was 
to calculate the equivalent yield-coupling factor. This was done by plotting the 
total energy in the Macor block as a function of time for each height of burst and 
depth of burial calculation. The total energy transmitted into the Macor block for 
the height of burst simulations is straightforward. The total energy in the Macor 
varies inversely with the height of burst. At the theoretical limit of an infinite 
height of burst, there would be no energy transferred into the Macor block, and 
the only energy in the Macor would be due to its inherent mass energy. Figure 27 
shows the energy transmitted into the Macor block as a function of time for each 
of the height of burst calculations. 
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Figure 27. Energy in Macor block for several HOB shots 


As expected, the energy decreases to a baseline as the height of burst 
increases. A similar relationship should also hold for the depth of burial shots. 
At the theoretical limit of an infinitely buried depth of burst, all the energy will 
remain within the Macor and none transmitted into the atmospheric gas mix. 
Figure 28 shows the energy contained within the Macor block as a function of 
time for the depth of burial calculations. 
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Figure 28. Energy in Macor block for several DOB shots 
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The total energy within the Macor increases from 5mm depth of burst to 
10mm depth of burst, yet there is a noticeable difference between the profiles of 
the 5mm and 10mm DOB energies and the others. As the depth of burial 
increases beyond 20mm DOB, the energy begins to decline with depth of burst. 
At first glance, it may seem to indicate an error in the data, as the energy should 
increase to an asymptote, as the depth of burial gets deeper. Upon further 
investigation, using LLNL’s visualization application, VisIT, it was discovered that 
for both the 5mm and 10mm depth of burial runs the exploding gas escapes the 
Macor block and leaves behind a crater. The point at which the Macor gives way 
corresponds to the two peaks in the total energy curve shown above for both the 
5mm and 10mm DOB shots, as shown in Figure 28. Figure 29 shows a split 
view movie frame, with material on the right and pressure on the left, at 35 
microseconds of the detonation. The breakout of the gold gas from the Macor 
block is clear. This breakout corresponds to the peak in the energy versus time 
plot in the Macor block for the 10mm DOB shot. 
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Figure 29. Split view of Material and Pressure at 35 microseconds for 10mm 

DOB shot 
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For all shots below 10mm, the fireball never escapes the Macor block, but 
simply bulges the “ground” where the fireball is forming. Figure 30 shows that 
the gold gas never escapes the Macor block by 120 microseconds for a 20mm 
depth of burial simulation. 



-J5# -IM 0 M J00 


Figure 30. Split view of Material and Pressure at 120 microseconds for 20mm 

DOB shot 

As demonstrated in Figures 30 and 31, the layer of Macor surrounding the 
gold gas becomes thicker as the depth of burial increases. 
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Figure 31. Split view of Material and Pressure at 120 microseconds for 30mm 

DOB shot 

Thus, for 20mm and lower depths of burial, the energy contained within 
the gold gas fireball was added to the total energy of the Macor block, and with 
this addition, the total energy asymptotes as expected. The baseline energy of 
the Macor, due to its mass energy, was subtracted from each to show that the 
energy asymptotes to the lOkJ value of total deposited energy. Figure 32 
illustrates the results of this corrected calculation of energy in the Macor block. 
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Figure 32. Corrected Plot of Energy in Macor block for several DOB shots 


With the corrected energy-time plots for both height of burst and depth of 
burial runs, the peak energy values were plotted as a function of depth of burial 
in Figure 33 in order to construct the more familiar energy coupling curve. 


Energy Coupling Cuive (Lab Frame) 



Figure 33. Lab Frame Energy Coupling Curve 
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In this plot, the ordinate corresponds to the energy coupling yield factor 
(ECYF), which is simply the peak energy in the Macor divided by the total input 
energy of 10kJ. A zero ECYF would correspond to an infinite high height of burst 
detonation and an ECYF of one would correspond to an infinitely deep depth of 
burial detonation. The abscissa indicates depth of burial where negative values 
along this axis correspond to the heights of burst for above-ground detonations. 
One key aspect of these data is that they are scalable. In order to scale up from 
the lab frame, in units of millimeters and kilojoules, the cube root scaling law, 
explained earlier in Chapter II sub-section B, “Scaling,” was used to convert the 
straight depth of burial to a scaled depth of burial with units of distance per yield 
energy^^^. The conversion factor shown below takes advantage of the cube root 
scaling law noted by Glasstone and Dolan. 

Conversion factor: 

Ifc/= 2.383*10(megaton) 

=(2.383*10=61.99*10 

1mm = 0.001m 

1mm / = 0.001m / (61.99 = 16.13m / MT^'^ 

Each of the depth of burial/ height of burst was scaled up to a scaled 
depth of burial/ scaled height of burst using the conversion factor. Flow each 
distance was scaled is shown below, using the 5mm depth of burial case as an 
example. 

simulation depth of burial/ yield'^^ * conversion factor = scaled depth of burial 

5mm/10fc/‘'^*(16.13m/Mr''^)/(lmm/lfc/‘'^) = 37.36m/Mr''^ 

Table 5 shows the results of similar calculations for each depth of burial 
ranging from 5mm to 50 mm with a 10kJ yield. 
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Table 5. 


rnmTTokJ^ 
5_ 

JO_ 

20 _ 

30_ 

40_ 

50 


_ m 

37.36 _ 

74.7 _ 

149.45 _ 

224.18 _ 

298.9 _ 

373.64 _ 

Scaled Depth of Burial 


Figure 34 shows the energy-coupling curve scaled to the world frame. 


Energy Coupling Curve Scaled Up to Real World Frame 



Figure 34. World Frame Energy Coupling Curve 

The data can now be mapped to earlier calculations of explosions in 
various hard rock media made by several researchers and summarized by L. 
Glenn at LLNL in Figure 35. Those calculations are focused in the regime where 
there is strong coupling between the ground and air shock. 
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Figure 35. Energy Coupling Curve for Hard Rock From L.GIenn, 2002, p.16 


First, the upper limit of the data fit was calculated by using the equation 
from the depth of burial data fit in Figure 35, setting the equivalent yield coupling 
factor (ECYF) to unity, and calculating the scaled depth of burial that 
corresponds to full coupling, i.e., no energy transmitted into air shock. 

1 = 0.22* scaled depth of burial (SDOB)° ^'’ 
ln(l/0.22) = 0.36*ln(5DOB) 

\niSDOB) = (1 / 0.36) * ln(l / 0.22) 

SDOB = 22) ^ 67_o8[m / 

This shows, according to the data fit, that all shots conducted at a scaled 
depth of burial greater than 67.08 m/ represent a situation in which all of 
the energy would be fully coupled into the ground. Using the scaled depth of 
burial points from Table 2 and the data fit formulas summarized by L. Glenn, the 

ECYF data points from the numerical simulation was plotted against the coupling 
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curve of the hard rock data in order to show the correlation with the Macor 
simulation model. Figure 36 shows that within the limits of the data fit model, 
both the simulation and the hard rock model data fit agree within a margin of 5%. 
The red triangular points are from the simulation data for the various HOB/DOB 
calculations. The green curve is a plot of the two formulas from the hard rock 
data equations. 


Matched Energy Coupling Curve 



Figure 36. Matched Coupling Curve 
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Table 6 shows the EYCF for each of the 12 simulation calculations. The 
EYCF was calculated by taking the peak energies, from Figures 27 and 32, for 
the corresponding height of burst and depth of burial calculations respectively, 
and then dividing that value by the lOkJ input energy. So, for example, in the 
10mm DOB calculation the maximum energy contained within the Macor was 
9.4kJ which corresponds to the peak of the 10mm curve in Figure 32. This value 
of 9.4kJ was then divided by the 10kJ yield to result in an ECYF of 0.94. 


Depth of Burial/ Fleight of Burst 

Equivalent Yield Coupling Factor 

5mm 

.85 

10mm 

.94 

20mm 

.99 

30mm 

.99 

40mm 

.99 

50mm 

.99 

-5mm (Height of Burst) 

.012 

-10mm 

.006 

-20mm 

.003 

-30mm 

.0025 

-40mm 

.002 

-50mm 

.0018 


Table 6. Equivalent Yield Coupling Factor from simulation data 


Table 7 shows a similar data chart, the difference being that the scaled 
depths of burial/ heights of burst (shown in Table 5) were used along with the two 
equations from the hard rock data for both the height of burst and depth of burial 
to calculate the EYCF. 
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Scaled DOB/ Scaled HOB [m/MT''^] 

Equivalent Yield Coupling Factor 

37.36 

.81 

74.7 

SDOB larger than 67.08 ^1 

149.45 

SDOB larger than 67.08 ^1 

224.18 

SDOB larger than 67.08 ^1 

298.9 

SDOB larger than 67.08 ^1 

373.64 

SDOB larger than 67.08 ^1 

-37.36 (Height of Burst) 

.012 (**outside data range yet agrees) 

-74.7 

.011 (outside data range) 

-149.45 

.010 (outside data range) 

-224.18 

.010 (outside data range) 

-298.9 

.010 (outside data range) 

-373.64 

.009 (outside data range) 


Table 7. Equivalent Yield Coupling Factor from live data fit curves 


The range in which both the simulation data and the hard rock data 
overlap is small, yet within that range both the simulation and hard rock data 
agree to within 5%. Even slightly beyond the limits of the curve fit from the hard 
rock data both the simulation and it agree, as noted in the 5mm HOB (-37.36 
m/MT^^^ scaled HOB), where both the simulation and the hard rock data fit predict 
the exact same ECYF of .012. 
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VII. CONCLUSION 


Scaled simulation and computer modeling of nuclear weapon effects will 
play a significant role in the maintenance of the SBSSP. The National Ignition 
Facility at LLNL provides a unique opportunity to investigate nuclear weapon 
effects and the physics of nuclear weapons in support of the achievement of 
goals put forth by the NNSA and the SBSSP. This work provides information in 
evaluating the hydrodynamic response of materials to nuclear detonations at or 
near the earth’s surface, which is a major component of the SPSSP agenda. 
The work presented in this thesis serves two purposes: first, to assist in the 
planning of the placement of shock sensors within a Macor object that will be 
used as a platform for upcoming NIF experiments designed to study the energy 
coupling/partitioning of a low yield nuclear detonation at or near the earth’s 
surface; and second, to develop a basic set of computational data for 
development of a coupling curve. The coupling curve derived here, though partly 
based on an idealistic buried device, has proven to be very consistent with 
calculations of previous explosions in various hard rock media supported by 
limited test data. The coupling curve derived here and the supporting 
calculations serve as a basic set of data for further investigation of buried nuclear 
devices. Using the work done in this thesis, opportunities for further research 
include adding the effect of radiation to the simulations and altering the design of 
the Macor object to more closely model an inhomogeneous material. 
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APPENDIX A 


Appendix A is the input file for a 10mm height of burst calculation. Many 
of the original marker files have been removed for the purpose of saving paper. 
This should help clarify particular portions of the text, specifically the topics of 
boundary conditions and equations of state. This input file was called into the 
GEODYN executable and the results analyzed with both MATLAB and VisIT. 


% units are specified in millimeters, milligrams, microseconds (GPa, km/s) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%% DEFINITIONS %%%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

HOB=10 % height of burst 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%% GLOBAL PARAMETERS %% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

max_step = 100000 
stop_time = 200 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%% DOMAIN GEOMETRY %%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

geometry:( 

% cylindrical (r-z) coordinates 
coord_sys = 1 

% low boundaries of (r,z) in domain 
probjo = [0 -120] 

% high boundaries of (r,z) in domain 
prob_hi = [160 160] 

) 


o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%%% MESH/OUTPUT %%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 


amr:( 

% number of elements in (r,z) on base grid 
n_cell = [80 140] 
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% verbose output 
V = 1 

% verbose timestep output 
v_step = 1 

% maximum level of refinement 
maxjevel = 2 

% regrid interval (timesteps) 
regridjnt = 2 
% maximum grid size 
max_grid_size = 32 
% refinement ratio (per level) 
ref_ratio = [44444] 

% number of extra cells around a grid to refine (per level) 
n_error_buf = [12111] 

% root name of restart files 
check_file = sesAucheck 
% restart file interval (steps) 
checkjnt = 100 

% root name of plot files 
plot_file = sesAuPlot 

% plot file interval (steps) [cannot be specified with plot_per] 

%plotJnt = 100 

% plot file period (microseconds) [cannot be specified with plotjnt] 
plot_per = 1 

% plot variables (see RULESETS) 

plot_vars = [ pressure density color vr vz sound_speed internal_energy 
bulking_porosity pressure_cutoff friction_slope T11 T22 T33 T12 temperature ] 

% marker specifications [ name (l)agrange/(e)ulerian r z ] 

Markers = [ 

[ pgaugel e 86.0 50.0 ] 

[ pgauge2 e 130.0 75.0] 

[ pgauge3 e 140.0 25.0] 

[gpgauge2 I 0.5 -70.0] 

[ gpgaugel I 35.0 -35.0 ] 

[gpgauge3 I 70.0 -70.0] 


] 

% marker output interval (steps) 
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markerjnt = 1 

% marker file write interval (steps) 
marker_dumpjnt = 100 
% marker variables (see RULESETS) 

marker_vars = [ m_pressure m_density m_vel_mag m_vr m_vz m_T11 m_T22 
m_T33 m_T12 m_color m_temperature mjnternal_energy m_soundspeed 
m_plastic_strain m_Tmin m_rdisp m_zdisp ] 

) 


o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%% PHYSICS/SOLVER %%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

hyp:( 

% number of materials 
n_materials=5 
% verbosity 

V = 0 

% minimum timestep below which to halt 
dt_cutoff = 0.00001 
% Courant stability factor 
cfl = .3 

% initial scale factor for timestep 
init_shrink = .01 

% maximum timestep increase per cycle 
change_max = 1.1 
% gravitational acceleration 
gravity = 0. 

% interval to sum integrated quantities 
sumjnterval = 100 

% integrated quantity variables (see RULESETS) 

integrated_quantities = [ sum_mass_source sum_energy_source 

sum_mass_source2 sum_energy_source2 sum_mass_source0 

sum_energy_source0 ] 

% write integrated quantities file on restart 
%i n teg rate_o n_restart=1 

) 


o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% PROBLEM CONDITIONS %% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

Prob:( 

% materials (see MATERIALS) 
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materials=[ gas_mix_LEOS Au_LEOS Macor Air_LEOS Fe_LEOS ] 

% boundary condition types 

%% 0 = Interior/Periodic 2 = Outflow %% 

%% 1 = Inflow 3 = Symmetry %% 

bc:lo:type = [31] 
bc:hi:type = [11] 

% inflow boundary specifications 
%% material is 0-indexed, density in g/cc, energy in kJ/g 
bc:lo:1:( 
material=3 
rho=0.0001 
e=0.52 

) 

bc:hi:0:( 

material=3 

rho=0.0001 

e=0.52 

) 

bc:hi:1:( 
material=3 
rho=0.0001 
e=0.52 

) 


% Specifies that displacement variables should be stored 
celldata<CellData>:( 
ops:[ 

<CellDisplacement>=0 

] 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%% CONFIGURATION %%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

Config:( 

% List of region names 

%% Specified order that regions are painted in. 

%% Last region fills remaining space. 

Regions = [ source air_cone steel_cone ceramic 
ring atmosphere ceramic_shell air ] 

% Region specifications 

%% material is 0-indexed, density in g/cc, energy in kJ/g 
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Region:( 

% source 

%% 10kJ in 2mg of material, sphere of diameter 2mm 
source:( 

% expanded iron 
material = 1 
rho = 0.4774648 
e = 5e+3 
file = circle512 

% no scaling, translate in +y by height of burst 
Transform = [ scale 1.1. translate 0. '$:HOB' ] 

) 

% air inside cone 
air_cone:( 

% air 

material = 3 
rho = 0.0001 
e = 0.52 

% r-z pairs specifying polygon (r-z projection of cone) 
Points = [ 

[0 '$:HOB+1000'] 

[ '0.6+1000/6' '$:HOB+1001'] 

[0.6'$:HOB+1' ] 

[0.'$:HOB' ] 

] 


% steel cone 
steel_cone;( 

% iron at reference density and energy 
material = 4 

% r-z pairs specifying polygon (r-z projection of outer cone) 
Points = [ 

['0.4+1000/6' '$:HOB+1001'] 

[0.4'$:HOB+1' ] 

[0.6'$:HOB+1' ] 

[ 1.4'$:HOB+1.8' ] 

['1.4+1000/6' '$:HOB+1001.8'] 

] 


% macor 
ceramic:( 

% macor at reference density 
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material = 2 
e = 0.0000405 
Points = [ 

[0 -100] 
[100 -100] 
[100 0.06] 
[0 0.06] 

] 


% macor ring 
ring:( 

% macor at reference density 
material = 2 
e = 0.0000405 


Points = [ 


[ 100.5 

-10] 

[ 149.5 

-10] 

[ 150. 

0.06] 

[ 100. 

0.06] 


] 


% gas-mixture atmosphere 
atmosphere:( 

% gas mixture 
material = 0 
rho=0.00153 
temperature = 310 
% composition of mixture 
"qq[0]" = 0.65 
"qq[1]" = 0.20 
"qq[2]" = 0.10 
"qq[3]" = 0.05 
Points = [ 

[0 - 100 ] 

[150 -100] 

[150 150] 

[ 0 150 ] 

] 

) 

% macor shell 
ceramic_shell:( 

% macor at reference density 
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material = 2 
e = 0.0000405 
Points = [ 

[0 -105] 

[155 -105] 

[155 155] 

[0 155] 

[0 150] 

[150 150] 

[150 -100] 

[0 - 100 ] 

] 

) 

% remaining air 
air:( 

% air 

material = 3 
rho = 0.0001 
e = 0.52 

) 

) 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%% REFINEMENT RULES %% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Specifies refinement condition 
refinement<Tagger>:( 

% refine on material interfaces and pressure gradients 
logic<SymbolFunction>:( 
f="(interfaces|pressure)" 
variable="interfaces,pressure" 

) 

% refinement variables (see RULESETS) 
derive=[ref_color ref_p ] 

% material interface detection 
%% refine when "color" is not an integer 
interfaces<FabFunction>:( 
comp=[ref_color] 

f=abs_remainder_offset parameters=[1 .e-8] 

) 
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% pressure gradient detection 

%% refine when relative difference (pA - pB)/(min(pA,pB,pMin)) > 0.3 
%% where pMin = 0.0002GPa 
pressure<FabFunction>:( 
comp=[ref_p] 

f=relative_gradient_offset parameters=[0.0002 .3] 

) 


) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%%%% RULESETS %%%%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% plot variables (rules must match names) 
pressure<DeriveRuleSet>:( 
models=[ StressTensor Energies ] 
rules:[ 

<Rule>:(var=p averager=bulk_modulus_average) 
<Rule>:(var=rho averager=volume_average) 
<Rule>:(var=material_number averager=volume_average) 
<Rule>:(var="v[0]") 

<Rule>:(var="v[1]") 

<Rule>:(var=c) 

<Rule>:(var=e) 

<Rule>:(var=bulking_porosity) 
<Rule>:(var=pressure_cutoff) 

<Rule>:(var=friction_slope) 

<Rule>:(var="T[T11]") 

<Rule>:(var="T[T22]") 

<Rule>:(var="T[T33]") 

<Rule>:(var="T[T12]") 

<Rule>:(var=temperature) 

] 

) 

% marker variables (rules must match names) 
m_pressure<DeriveRuleSet>:( 
models=[StressTensor Energies ] 
rules:[ 

<Rule>:(var=p averager=bulk_modulus_average) 
<Rule>:(var=rho averager=volume_average) 
<Rule>:(var=v_mag) 

<Rule>:(var="v[0]") 

<Rule>:(var="v[1]") 
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<Rule>:(var="T[T11]") 

<Rule>:(var="T[T22]") 

<Rule>:(var="T[T33]") 

<Rule>:(var="T[T12]") 

<Rule>:(var=material_number averager=volume_average) 
<Rule>:(var=temperature) 

<Rule>:(var=e) 

<Rule>:(var=c) 

<Rule>:(var=plastic_strain) 

<Rule>:(var=Tminhistory) 

<Rule>:(var="displacement[0]") 
<Rule>:(var="displacement[1 ]") 

] 


% integrated quantities variables (rules must match names) 
sum_mass_source<DeriveRuleSet>:( 
models = [ Energies ] 
rules:[ 

<Rule>:(var=rho averager=volume_sum materials=[1]) 
<Rule>:(var=e_tot averager=mass_sum materials=[1]) 
<Rule>:(var=rho averager=volume_sum materials=[2]) 
<Rule>:(var=e_tot averager=mass_sum materials=[2]) 
<Rule>:(var=rho averager=volume_sum materials=[0]) 
<Rule>:(var=e_tot averager=mass_sum materials=[0]) 

] 

) 

% refinement variables 
ref_color<DeriveRuleSet>:( 
models=[EOS] 
rules;[ 

<Rule>:(var=material_number averager=volume_average) 
<Rule>:(var=p averager=bulk_modulus_average) 

] 

) 


o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%%%% MATERIALS %%%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

gas_mix_LEOS<Material>: ( 
models: [ 

<NoPorosityVariable> = 0 
<RhoSolid> = 0 
<N_GammaAveraged_EOS>:( 
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EOS:[ 

<LEOS> :( use_energy_offset=1 material_number=100 ) % qq[0] - Ne 
<LEOS> ;( use_energy_offset=1 material_number=180 ) % qq[1] - Ar 
<LEOS> :( use_energy_offset=1 material_number=360 ) % qq[2] - Kr 
<LEOS> :( use_energy_offset=1 material_number=540 ) % qq[3] - Xe 

] 

) 

<N_LimitFraction> = 0 
<SphericalStressTensor> = 0 
<ZeroShearModulus> = 0 
<Smooth_PressureCutOff>: (pmin=0) 

<ZeroFrictionSlope> = 0 
<Tminhistory> = 0 

] 

) 

Air_LEOS<Material>: ( 
models: [ 

<NoPorosityVariable> = 0 
<RhoSolid> = 0 

<LEOS> :( material_number=2260 ) 

<NoUpdate> = 0 
<SphericalStressTensor> = 0 
<ZeroShearModulus> = 0 
<Smooth_PressureCutOff>: (pmin=0) 

<ZeroFrictionSlope> = 0 
<Tminhistory>=0 

] 

) 

Fe_LEOS<Material>: ( 
models: [ 

<NoPorosityVariable> = 0 
<RhoSolid> = 0 

<LEOS>: ( material_number=260 ) 

<NoUpdate> = 0 
<SphericalStressTensor> = 0 
<ZeroShearModulus> = 0 
<Smooth_PressureCutOff>: (pmin=0) 

<ZeroFrictionSlope> = 0 
<Tminhistory>=0 

] 

) 

Au_LEOS<Material>: ( 
models: [ 
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<NoPorosityVariable> = 0 
<RhoSolid> = 0 

<LEOS>: ( material_number=2700 leos_name=/usr/gapps/data/eos/sesame ) 

<NoUpdate> = 0 

<SphericalStressTensor> = 0 

<ZeroShearModulus> = 0 

<Smooth_PressureCutOff>: (pmin=0) 

<ZeroFrictionSlope> = 0 
<Tminhistory>=0 

] 

) 


Macor<Material>: ( 
models: [ 

<ElasticPorosityVariable> = 0 
<RhoSolid> = 0 
<MieGruneisen_EOS>: ( 
rhoO = 2.52 
cO =4.58 
s1 =1.4 
s2 =0. 

S3 =0. 

g = 1 - 

b = 0. 

) 

<HistoryDependent_PressureCutOff>: ( pmin=-0.008) 
<PoreElasticity>: ( a=0. b=1.) 

<PseudoCapYieldStrength1 >:( 
scale_factor=1. 

Y_c=0.345 

Y_t=0.1 

tau_dam=0.00001 
P_hard_rate=5. 

P_hard_exp=0.5 

Yrat=0.98 

cap0=0.5 

cap00=0.5 

cap_power=2. 

strain_hard_eps=0.02 

strain_tofail=0.0001 

soft_rate=5. 

rate_function<SymbolicFunction>: (function="(x+1 .)'^0.12") 
residual_function<SymbolicFunction>: (function="1 .*x") 

) 
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<PseudoCapCompaction>:( 

poro0=0.001 

mu_c0=10. 

slope0=0.1 

bulking_k1<SymbolicFunction>: (function="0.05*x") 
tau_comp=0 

) 

<RateDependentFunctionPlasticUpdatewithCorrection>:( 
tau_of_p<Const>: (f=0 ) 

) 

<ElasticUpdate1> = 0 

<YieldLodeAngleTerm>: ( q20 = 0.5 function<Const>:(f=1)) 

<Stress> = 0 

<ConstantPoissonRatio>: ( nu=0.29 ) 

<PseudoCapBulkingwithCorrection> = 0 

<PseudoCapBulkingCorrection>: ( 
k1=0.5 
k=0.5 

) 

<FrictionSlope>: ( y0=1.5 ) 

] 
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APPENDIX B 


Appendix B is the input file for the modified 30mm depth of burial 
calculation. Many of the original marker files have been removed for the purpose 
of saving paper. The original goal of the modified DOB set up was to have a 
layer of air around both the Macor and the gas-mix, but without a solid boundary 
between the two gases the gas-mix expanded to fill the problem space. This 
setup worked also as it still left a single material at the problem boundary, which 
eliminated the reflected wave noise, as shown in Figures 17 and 18. This input 
file was called into the GEODYN executable and the results analyzed with both 
MATLAB and VisIT. 

% units are specified in millimeters, milligrams, microseconds (GPa, km/s) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%% DEFINITIONS %%%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

FIOB=-30 % depth of burst 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%% GLOBAL PARAMETERS %% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

max_step = 100000 
stop_time = 200 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%% DOMAIN GEOMETRY %%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

geometry:( 

% cylindrical (r-z) coordinates 
coord_sys = 1 

% low boundaries of (r,z) in domain 
probjo = [0 -120] 

% high boundaries of (r,z) in domain 
prob_hi = [160 160] 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%%% MESH/OUTPUT %%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

amr:( 
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% number of elements in (r,z) on base grid 
n_cell = [80 140] 

% verbose output 
V = 1 

% verbose timestep output 
v_step = 1 

% maximum level of refinement 
maxjevel = 2 

% regrid interval (timesteps) 
regridjnt = 2 
% maximum grid size 
max_grid_size = 32 
% refinement ratio (per level) 
ref_ratio = [44444] 

% number of extra cells around a grid to refine (per level) 
n_error_buf = [12111] 

% root name of restart files 
check_file = chk 
% restart file interval (steps) 
check int = 200 


% root name of plot files 
plot_file = pit 

% plot file interval (steps) [cannot be specified with plot_per] 

%plotJnt = 100 

% plot file period (microseconds) [cannot be specified with plotjnt] 
plot_per = 1 

% plot variables (see RULESETS) 

plot_vars = [ pressure density color vr vz sound_speed internal_energy 
bulking_porosity pressure_cutoff friction_slope T11 T22 T33 T12 temperature ] 


% marker specifications [ name (l)agrange/(e)ulerian r z ] 
Markers = [ 

86.0 50.0 ] 

130.0 75.0] 

25.0 140.0] 

0.5 -70.0 ] 

I 35.0 -35.0 ] 

I 70.0 -70.0 ] 


[pgaugel 
[pgauge2 
[pgauge3 
[gpgauge2 
[gpgaugel 
[gpgauge3 


e 

e 

e 

I 
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] 

% marker output interval (steps) 
markerjnt = 1 

% marker file write interval (steps) 
marker_dumpjnt = 100 
% marker variables (see RULESETS) 

marker_vars = [ m_pressure m_density m_vel_mag m_vr m_vz m_T11 m_T22 
m_T33 m_T12 m_color m_temperature mjnternal_energy m_soundspeed 
m_plastic_strain m_Tmin m_rdisp m_zdisp ] 

) 


o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%% PHYSICS/SOLVER %%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

hyp:( 

% number of materials 
n_materials=4 
% verbosity 

V = 0 

% minimum timestep below which to halt 
dt_cutoff = 0.00001 
% Courant stability factor 
cfl = .2 

% initial scale factor for timestep 
init_shrink = .01 

% maximum timestep increase per cycle 
change_max = 1.1 
% gravitational acceleration 
gravity = 0. 

% interval to sum integrated quantities 
sumjnterval = 100 

% integrated quantity variables (see RULESETS) 

integrated_quantities = [ sum_mass_source sum_energy_source 

sum_mass_source2 sum_energy_source2 sum_mass_source0 

sum_energy_source0 ] 

% write integrated quantities file on restart 
%i n teg rate_o n_restart= 1 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% PROBLEM CONDITIONS %% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

Prob:( 
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% materials (see MATERIALS) 

materials=[ gas_mix_LEOS Au_LEOS Macor Air_LEOS ] 

% boundary condition types 

%% 0 = Interior/Periodic 2 = Outflow %% 

%% 1 = Inflow 3 = Symmetry %% 

bc:lo:type = [31] 
bc:hi:type = [11] 

% inflow boundary specifications 
%% material is 0-indexed, density in g/cc, energy in kJ/g 
bc:lo:1:( 
material=3 
rho=0.0001 
e=0.52 

) 

bc:hi:0:( 

material=3 

rho=0.0001 

e=0.52 

) 

bc:hi:1:( 
material=3 
rho=0.0001 
e=0.52 

) 


% Specifies that displacement variables should be stored 
celldata<CellData>:( 
ops:[ 

<CellDisplacement>=0 

] 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%% CONFIGURATION %%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

Config:( 

% List of region names 

%% Specified order that regions are painted in. 

%% Last region fills remaining space. 

Regions = [ source ceramic atmosphere ] 

% Region specifications 
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%% material is 0-indexed, density in g/cc, energy in kJ/g 
Region:( 

% source 

%% 10kJ in 2mg of material, sphere of diameter 2mm 
source:( 

% expanded gold 
material = 1 
rho = 0.4774648 
e = 5e+3 
file = circle512 

% no scaling, translate in +y by height of burst 
Transform = [ scale 1.1. translate 0. '$:HOB' ] 

) 

% macor 
ceramic:( 

% macor at reference density 
material = 2 
e = 0.0000405 
Points = [ 

[0 - 110 ] 

[150 -110] 

[150 0.06] 

[0 0.06] 

] 

) 

% gas-mixture atmosphere 
atmosphere:( 

% gas mixture 
material = 0 
rho=0.00153 
temperature = 310 
% composition of mixture 
"qq[0]" = 0.65 
"qq[1]" = 0.20 
"qq[2]" = 0.10 
"qq[3]" = 0.05 
Points = [ 

[0 0.06] 

[0 150] 

[150150] 

[150 0.06] 

] 

) 
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% air boundary 
air:( 

% air boundary 
material = 3 
e = 0.52 
rho=0.0001 

) 


) 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%% REFINEMENT RULES %% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Specifies refinement condition 
refinement<Tagger>:( 

% refine on material interfaces and pressure gradients 
logic<SymbolFunction>:( 
f="(interfaces|pressure)" 
variable="interfaces,pressure" 

) 

% refinement variables (see RULESETS) 
derive=[ref_color ref_p ] 

% material interface detection 
%% refine when "color" is not an integer 
interfaces<FabFunction>:( 
comp=[ref_color] 

f=abs_remainder_offset parameters=[1 .e-8] 

) 

% pressure gradient detection 

%% refine when relative difference (pA - pB)/(min(pA,pB,pMin)) > 0.3 
%% where pMin = 0.0002GPa 
pressure<FabFunction>:( 
comp=[ref_p] 

f=relative_gradient_offset parameters=[0.0002 .3] 

) 


) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

O/ O/ O/ O/ O/ O/ RIM RQRTQ O/ O/ O/ O/ O/ O/ O/ 

/o /o /o /o /o /o nULCOC I O /o /o /o /o /o /o /o 
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o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 


% plot variables (rules must match names) 
pressure<DeriveRuleSet>:( 
models=[ StressTensor Energies ] 
rules:[ 

<Rule>:(var=p averager=bulk_modulus_average) 
<Rule>:(var=rho averager=volume_average) 
<Rule>:(var=material_number averager=volume_average) 
<Rule>:(var="v[0]") 

<Rule>:(var="v[1]") 

<Rule>:(var=c) 

<Rule>:(var=e) 

<Rule>:(var=bulking_porosity) 
<Rule>:(var=pressure_cutoff) 

<Rule>:(var=friction_slope) 

<Rule>:(var="T[T11]") 

<Rule>:(var="T[T22]") 

<Rule>:(var="T[T33]") 

<Rule>:(var="T[T12]") 

<Rule>:(var=temperature) 

] 

) 

% marker variables (rules must match names) 
m_pressure<DeriveRuleSet>:( 
models=[StressTensor Energies ] 
rules:[ 

<Rule>:(var=p averager=bulk_modulus_average) 
<Rule>:(var=rho averager=volume_average) 
<Rule>:(var=v_mag) 

<Rule>:(var="v[0]") 

<Rule>:(var="v[1]") 

<Rule>:(var="T[T11]") 

<Rule>:(var="T[T22]") 

<Rule>:(var="T[T33]") 

<Rule>:(var="T[T12]") 

<Rule>:(var=material_number averager=volume_average) 
<Rule>:(var=temperature) 

<Rule>:(var=e) 

<Rule>:(var=c) 

<Rule>:(var=plastic_strain) 

<Rule>:(var=Tminhistory) 

<Rule>:(var="displacement[0]") 
<Rule>:(var="displacement[1 ]") 
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] 

) 

% integrated quantities variables (rules must match names) 
sum_mass_source<DeriveRuleSet>:( 
models = [ Energies ] 
rules:[ 

<Rule>:(var=rho averager=volume_sum materials=[1]) 
<Rule>:(var=e_tot averager=mass_sum materials=[1]) 
<Rule>:(var=rho averager=volume_sum materials=[2]) 
<Rule>:(var=e_tot averager=mass_sum materials=[2]) 
<Rule>:(var=rho averager=volume_sum materials=[0]) 
<Rule>:(var=e_tot averager=mass_sum materials=[0]) 


] 

) 

% refinement variables 
ref_color<DeriveRuleSet>:( 
models=[EOS] 
rules:[ 

<Rule>:(var=material_number averager=volume_average) 
<Rule>:(var=p averager=bulk_modulus_average) 

] 

) 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

%%%%%% MATERIALS %%%%%% 

o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ o/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

gas_mix_LEOS<Material>: ( 
models: [ 

<NoPorosityVariable> = 0 
<RhoSolid> = 0 
<N_GammaAveraged_EOS>:( 

EOS:[ 

<LEOS> :( use_energy_offset=1 material_number=100 ) % qq[0] - Ne 
<LEOS> :( use_energy_offset=1 material_number=180 ) % qq[1] - Ar 
<LEOS> :( use_energy_offset=1 material_number=360 ) % qq[2] - Kr 
<LEOS> :( use_energy_offset=1 material_number=540 ) % qq[3] - Xe 

] 

) 

<N_LimitFraction> = 0 
<SphericalStressTensor> = 0 
<ZeroShearModulus> = 0 
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<Smooth_PressureCutOff>: (pmin=0) 
<ZeroFrictionSlope> = 0 
<Tminhistory> = 0 

] 


Au_LEOS<Material>: ( 
models: [ 

<NoPorosityVariable> = 0 
<RhoSolid> = 0 

<LEOS>: ( material_number=2700 leos_name=/usr/gapps/data/eos/sesame ) 

<NoUpdate> = 0 

<SphericalStressTensor> = 0 

<ZeroShearModulus> = 0 

<Smooth_PressureCutOff>: (pmin=0) 

<ZeroFrictionSlope> = 0 
<Tminhistory>=0 

] 

) 

Air_LEOS<Material>:( 
models: [ 

<NoPorosityVariable>=0 

<RhoSolid>=0 

<LEOS>:(material_number=2260) 

<NoUpdate>=0 
<SphericalStressTensor>=0 
<ZeroShearModulus>=0 
<Smooth_PressureCutOff>: (pmin=0) 

<ZeroFrictionSlope>=0 

<Tminhistory>=0 

] 

) 


Macor<Material>: ( 
models: [ 

<ElasticPorosityVariable> = 0 
<RhoSolid> = 0 
<MieGruneisen_EOS>: ( 
rhoO = 2.52 
cO =4.58 
s1 =1.4 
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) 


s2 =0. 
S3 =0. 


g = 1 - 

b = 0. 

) 

<HistoryDependent_PressureCutOff>: ( pmin=-0.008) 
<PoreElasticity>: ( a=0. b=1.) 

<PseudoCapYieldStrength1 >:( 
scale_factor=1. 

Y_c=0.345 

Y_t=0.1 

tau_dam=0.00001 
P_hard_rate=5. 

P_hard_exp=0.5 

Yrat=0.98 

cap0=0.5 

cap00=0.5 

cap_power=2. 

strain_hard_eps=0.02 

strain_tofail=0.0001 

soft_rate=5. 

rate_function<SymbolicFunction>: (function="(x+1 .)''0.12") 
residual_function<SymbolicFunction>: (function="1 .*x") 

) 

<PseudoCapCompaction>:( 

poro0=0.001 

mu_c0=10. 

slope0=0.1 

bulking_k1<SymbolicFunction>: (function="0.05*x") 
tau_comp=0 

) 

<RateDependentFunctionPlasticUpdatewithCorrection>:( 
tau_of_p<Const>: (f=0) 

) 

<ElasticUpdate1> = 0 

<YieldLodeAngleTerm>: ( q20 = 0.5 function<Const>:(f=1)) 
<Stress> = 0 

<ConstantPoissonRatio>: (nu=0.29) 
<PseudoCapBulkingwithCorrection> = 0 
<PseudoCapBulkingCorrection>: ( 
k1=0.5 
k=0.5 

) 

<FrictionSlope>: ( y0=1.5 ) 
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